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ABSTRACT 

We use numerical simulation of cosmological reionization to calculate the secondary CMB 
anisotropies in a representative flat low density cosmological model. We show that the kinetic 
Sunyaev-Zel'dovich effect (scattering off of moving electrons in the ionized intergalactic medium) 
is dominated by the nonlinear hydrodynamic and gravitational evolution of the density and 
velocity fields, rather than the detailed distribution of the ionization fraction ("patchy reion- 
ization") on all angular scales. Combining our results with the recent calculation of secondary 
CMB anisotropies by Springel et al., we are able to accurately predict the power spectrum of the 
kinetic SZ effect on almost all angular scales. 

Subject headings: cosmic microwave background - cosmology: theory - cosmology: large-scale structure 
of universe - galaxies: formation - galaxies: intergalactic medium 



1. Introduction 

The reionization of the intergalactic medium 
(hereafter, IGM) stands out as one of the most im- 
portant physical processes that have taken place in 
the early universe. Recent observational progress 
measuring the angular spectrum of anisotropies 
in the Cosmic Microwave Background radiation 
(hereafter CMB) on sub-degree angular scales (de 
Bernardis et al. 2000; Hanany et al. 2000), and 
in obtaining the absorption spectra of the most 
distant quasars (Stern et al. 2000; Zheng et al. 
2000; Fan et al. 2000a) allows us to limit the red- 
shift of reionization to somewhere between 6 and 
30 (Griffiths, Barbosa, & Liddle 1998; Bond & 
Jaffe 1999; Tegmark & Zaldarriaga 2000). Com- 
bined with recent theoretical and numerical break- 
throughs in our understanding how reionization 
proceeds in the inhomogeneous universe (Giroux 
& Shapiro 1996; Haiman & Loeb 1997; Gnedin 
& Ostriker 1997; Madau, Haardt, & Rees 1999; 
Miralda-Escudee, Haehnelt, & Rees 1999; Ciardi 
et al. 1999; Valageas & Silk 1999; Chiu & Os- 



triker 2000; Gnedin 2000), these data finally place 
a study of reionization on much more solid footing. 

Among the many effects of reionization, the 
secondary CMB anisotropies have been the sub- 
ject of active study for a long time, beginning 
with the pioneering works of Ostriker & Vish- 
niac (1986) and Vishniac (1987) and up to the 
most recent careful investigations (Hu & White 
1996; Aghanim et al. 1996; Knox, Scoccimarro, 
& Dodelson 1998; Gruzinov & Hu 1998; Jaffe & 
Kamionkowski 1998; Peebles & Juszkiewicz 1998; 
Haiman & Knox 1999; Hu 2000; Bruscoh et al. 
2000; da Silva et al. 2000; Refregier et al. 2000; 
Seljak, BurweU, & Pen 2000; Springel, White, & 
Hernquist 2000). However, most of the recent 
studies have not benefited from the recent progress 
on reionization, and were based on over-simplified 
ad hoc models. In addition, many of these studies 
only considered the formerly favored "Standard" 
Cold Dark Matter (CDM) model. However, there 
appears finally to be convergence toward a low 
density CDM model with the cosmological con- 
stant (the so-called CDM-I-A model; c.f. White, 



Scott, & Picrpaoli 2000; Tegmark & Zaldarriaga 
2000; Bridle et al. 2000; Hu et al. 2000; Jaffe et al. 
2000), and thus it would make sense to calculate 
the secondary anisotropies for this currently most 
favored model. 

In this paper we reconsider the effect of reion- 
ization on the CMB, using the latest numerical 
simulations of reionization (Gnedin 2000). By 
combining our numerical results with analytical 
calculations and with recent numerical calculation 
by Springel et al. (2000) on larger angular scales, 
we are able to extend the range of angular scales 
over which predictions are reliable to over more 
than six orders of magnitude. We also take spe- 
cial care to make sure that our results are not con- 
taminated by numerical artifacts, the dominant of 
which is the periodicity of the simulation bound- 
ary conditions. 

2. Theory 

2.1. Secondary CMB Anisotropies 

Thomson scattering of CMB photons by ionized 
gas is usually called the Sunyaev-Zel'dovich (SZ) 
effect, which is often divided into "thermal" and 
"kinetic" SZ effects. The "thermal" SZ effect is 
the upscattering of CMB photons after compton 
scattering by hot gas (most pronounced in clus- 
ters of galaxies). The "kinetic" SZ effect is the 
temperature fluctuation induced by bulk motions 
of the gas. Both effects play a role in generat- 
ing secondary CMB anisotropies during and after 
reionization. 

The fractional temperature perturbation in the 
direction 6 induced by bulk motions ("kinetic ef- 
fect") is 
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where rig is the electron density along the line of 
sight, v{w; w) is the bulk peculiar velocity at posi- 
tion at a conformal time rjo — w, ar is the Thom- 
son cross section, t is the optical depth from us to 
w, and the subscript denotes a quantity at the 
current moment. The factor of a(w) arises because 
the physical time is dt = a{w)dw. Note that rep- 
resents a three-dimensional unit vector along the 
line of sight, whereas 6 will refer to a dimensionless 
two-dimensional vector in the plane perpendicu- 



lar to the line of sight — i.e., for directions near z, 
e = [91,02, ^/l- O'i-0'i) - (6*1, 6*2,1), whereas 
e= {01,02,0). 

It is convenient to rewrite equation (1) in the 
following dimensionless form: 
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where tq = cr]ocrTinH,o + 2nHe,o) and 
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Here Xe is the ionization fraction, S is cosmic over- 
density, and fin and nne are the mean hydrogen 
and helium number densities. 

Analogously, the "thermal" SZ effect is de- 
scribed by the following integral (for observations 
in the Rayleigh- Jeans regime): 

AT r« _,M7V-7L) dw 

where Tg and Tq are the gas and the CMB tem- 
peratures respectively. 

For a particular realization, we can simply cal- 
culate integral (1) or (3) along each line of sight 
and make a map of temperature anisotropy. In 
this paper we use numerical simulation of cosmo- 
logical reionization (described in detail in the fol- 
lowing section) to directly calculate the anisotropy 
map on the patch on the sky subtended by the 
eoniputational box. 

Often however it is just the two-point cor- 
relation function of the anisotropies, C{9), or 
its spherical-harmonic transform, the power spec- 
trum, Ci, which arc of interest. In this case we 
can directly calculate the correlation function from 
equation (1), 

C{0) ^ {^{0,)^{0,))^^^^^g^,g^, (4) 

which can be reduced to the following form: 
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where Cee and Cw are the electron density and 
velocity correlation functions, 

Cee{vi,m,^) = {xiri,rji)x{r2,m)), 
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with cos 6 = (r*! • r2)/(rir2), provided wc make a 
usual assumption (c.f. Bruscoli et al. 1999) that 
the electron density and velocity are uncorrelated, 
i.e., 

Ce.ivuV2,0) = {x{rum)02 ■ lir2,m)) = 0. (7) 

We will return to the accuracy of this assumption 
below. 

2.2. Terminology for Secondary Anisotropies 

In reality we have to deal with only one uni- 
verse, and thus we can only observe the total 
CMB anisotropies given by equation (1), gcneri- 
cally labeled as the kinetic Sunyaev-Zel'dovich ef- 
fect. Historically, this was first investigated by 
restricting to a much simpler, unphysical situa- 
tion: homogeneous reionization and the linear evo- 
lution of velocity and density fields. This is the 
"Ostriker-Vishniac effect" (OV; Ostriker & Vish- 
niac 1986; Jaffe & Kamionkowski 1998), also some- 
times called the "Vishniac effect." In this case, we 
observe the autocorrelation of the total gas density 
and velocity fields, related to the power spectrum 
of density fluctuations by linear theory. 

In this work, we will further distinguish be- 
tween the "Linear Ostriker-Vishniac Effect" (LOV) 
and the "Nonlinear Ostriker-Vishniac effect" 
(NLOV). The former is the usual calculation of 
OV, using the linear relationship between velocity 
and density. This comes into the needed calcu- 
lation of the four-point fimction (vSvS), in which 
it is assumed that (1) the linear relationship be- 
tween V and 6 holds; and (2) that the four-point 
function is appropriate for an underlying Gaus- 
sian density field, which only holds in the fully 
linear regime (and only with Gaussian initial con- 
ditions). The NLOV assumes a homogeneous ion- 
ization fraction, but the full nonlinear gas density 
and velocity fields. 

Recently, the LOV calculations have been com- 
plimented by consideration of the distribution of 
the ionized hydrogen in the universe as it under- 
goes reionization, so-called "patchy" or "inhomo- 
geneous" reionization. Workers in the field have 
considered models for the distribution of the ion- 
ization fraction in the linear regime (e.g., Knox, 
Scoccimarro, & Dodelson 1998; Gruzinov & Hu 
1998), as well as full simulations or semianalytic 
models of the ionized gas, as in this work (see also 



Bruscoli et al. (2000), Benson et al. (2000), and 
Springel et al (2000)). 

This separation is obviously somewhat artifi- 
cial. Nevertheless, for academic purposes, we can 
denote those two effects in the following fashion: 
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for the NL Ostriker-Vishniac effect, and 
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for patchy reionization. 

Thus, the NL Ostriker-Vishniac effect describes 
anisotropies generated in the universe in which the 
ionization fraction is homogeneous in space ("ho- 
mogeneous reionization"), and patchy reionization 
includes the rest of the anisotropies. We again 
emphasize that this separation is artificial and un- 
physical, and we adopt it purely for historical rea- 
sons, as the two effects were often considered to 
be "distinct" effects. 

Using the above definition, we can also define 
the electron density correlation functions due to 
NL Ostriker-Vishniac and patchy reionization ef- 
fects, 

CP^iVuV2,0) ^ (Ax(fi,77i)Ax(r2,772)XlO) 

where Ax = (xe — .Te)(l+(5). The total correlation 
function Cee is then a sum of these two functions 
plus the cross-correlation. 
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3. Method 
3.1. Simulation 

We use a cosmological simulation of reioniza- 
tion reported in Gnedin (2000). The simulation 
includes 3D radiative transfer (in an approximate 
implementation) and other physical ingredients re- 
quired for modeling the process of cosmological 
reionization. 
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The simulation of a representative CDM+A 
cosmological model^ was performed in a comov- 
ing box with the size of 4/i~^ Mpc with the mass 
resolution of 5 x 10''' Mq in baryons and the co- 
moving spatial resolution of lh~^ kpc. 

The simulation was stopped at 2 = 4 because at 
this time the rms density fluctuation in the com- 
putational box is about 0.25, and at later times 
the box ceases to be a representative region of the 
universe. 

As was noted in Gnedin (2000), this simula- 
tion is still insufficient to give the fully converged 
numerically result. As we show below, the com- 
putational box size of this simulation is not large 
enough to allow for accurate computation of the 
CMB anisotropics, as perturbations at ^; < 4 
contribute about 20% of the signal. Wc there- 
fore must emphasize here that our calculations are 
valid only on a semi-qualitative basis, within a fac- 
tor of 1.5-2, and a still larger simulation is required 
to accurately model all the relevant scales present 
in the problem. But since the observations of the 
secondary anisotropics are some years away, theo- 
rists have time to improve upon their models and 
perform larger, highly accurate simulations. 

Our simulation assumes that all cosmologi- 
cal reionization occurs via radiation from stel- 
lar sources, with star formation parameterized by 
the phenomenological Schmidt law as discussed in 
Gnedin (2000). We ignore alternative and compli- 
mentary scenarios in which the bulk of reionization 
is produced by Active Galactic Nuclei. In light of 
recent high-redshift quasar counts (e.g.. Fan et 
al. 2000b), a scenario in which the universe is 
rcionizcd by bright optically selected QSOs seems 
unlikely in any event. There remains however the 
possibility that low brightness AGNs contributed 
significantly (if not dominantly) to the reioniza- 
tion of the universe. In the latter case however 
they will be clustered on scales similar to the 
stellar sources, and thus will not result in a qual- 
itatively different reionization scenario, although 
the epoch of reionization in this case is less well- 
determined. 



With the following cosmological parameters: f!o = 0.3, 
= 0.7, h = 0.7, Ut = 0.04, n = 1, ag = 0.91, where 
the amplitude and the slope of the primordial spectrum are 
fixed by the COBE and cluster-scale normalizations. 



3.2. Numerical Issues 

Given a cosmological simulation, wc arc now 
faced with the task of calculating the secondary 
CMB anisotropics accurately from the simulation. 
Since perturbations are generated over a consider- 
able redshift range, ideally we would like to have 
a simulation with the box size that extends from 
z = until z ^ 30 in one direction. However, with 
currently feasible simulations this is impossible, 
and we utilize the cubic box with the 4/i~^ Mpc 
side used in Gnedin (2000). 

If we take the output of a simulation as is, the 
artificial periodicity of the universe will amplify 
the anisotropics by a large factor, simply because 
in real universe the signal will be averaged over 
many randomly placed Hll regions, while in a sim- 
ulated periodic universe a photon, on its way along 
the line of sight, will encounter the same struc- 
ture over and over again. This case is depicted in 
the first two rows of Figure 1. In order to avoid 
this artificial amplification, we randomly flip and 
transpose the computational box around any of 
its six edges, and in addition shift it by a random 
distance in a random direction. This "random- 
ization" procedure ensures that the next periodic 
image of the computational box is not correlated 
with the given image, and thus we lose the corre- 
lations over the scales larger than half the size of 
the computational box along the line-of-sight di- 
rection. (In addition to the flipping, transposition 
and shifting of the simulation box, we have also 
rotated it by a random amount, but this has neg- 
ligible effect). This is sketched in the third row 
of Fig. 1. Thus, our procedure actually underes- 
timates true anisotropics by ignoring large-scale 
signal. However, we can add the missing large- 
scale signal using the linear theory calculation, as 
we discuss below. 

In addition, since we cannot unambiguously 
eliminate periodicity in the plane of the sky, we 
restrict our image to the image of one periodic 
simulation box, which corresponds to the angular 
size of 

6 - ^ 

where L = 4ft.~^ Mpc is the comoving size of our 
computational box, R{t) is the comoving angular- 
diameter distance from time t to the current mo- 
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Fig. 1. — Cartoon showing the effect of the pe- 
riodic box on the calculated CMB anisotropics. 
Right panels depict a piece of the universe along 
the line of sight with Hii regions in it, and left 
panels show the correlation function of any physi- 
cal quantity as a function of the distance w along 
the line of sight. The upper row depicts the ac- 
tual Universe, with random positions for Hii re- 
gions and the correlation function extending to 
large scales. The middle row shows the periodic 
universe as given in a cosmological simulation. Fi- 
nally, the lower row depicts our method of calcu- 
lating the anisotropics, with the periodic box ran- 
domly shifted and flipped in the subsequent im- 
ages. Arrows in the last two rows show the orien- 
tation of the computational box. For simplicity, 
we do not show random shift of the box in the last 
row which is actually performed in our calculation 
in addition to random flipping and transposing. 



ment, 

and we choose to scale all the computational boxes 
along the line-of-sight to the box taken at aim = 
a(iim) = l/(l + ^im) = 0.115, the epoch just before 
reionization which contributes most to the CMB 
anisotropics. For our calculations, ^max — 60 arc- 
sec corresponding to a spherical harmonic multi- 
pole of 10^. 

Because in the simulation the bulk velocity of 
the computational box is set to zero, whereas a 
volume with the size of our box in real universe will 
have some non-vanishing bulk velocity, we add to 
each box a random bulk velocity drawn from the 
Gaussian distribution with the dispersion equal to 
the linear rms velocity on the scale of the compu- 
tational box at each redshift. Because this effect 
depends on the correlation of density with velocity, 
our results on all angular scales depend on such an 
appropriate realization of the velocities. 

In summary, our procedure to calculate the 
CMB anisotropics from the cosmological simula- 
tion with a periodic computational box can be de- 
scribed as follows: 

1. We select the outputs from the simulation 
spaced with Aa in scale factor. 

2. For each pair of outputs at and a^+i, 
where index i numbers the outputs, we cal- 
culate the number of computational boxes 
along the line of sight between and a^+i. 

3. Then for all boxes between and (oj -|- 
aj+i)/2 we fill them in with the simulation 

output data at a^, and the boxes between 
(oj -|- ai+i)/2 and Oi+i we fill in with the 
data at Oj+i. 

4. Then each of thus-filled-in boxes (including 
the original simulation boxes at ai and Oi+i) 
is randomized by a random shift in the plane 
of the sky and by a random flip and transpo- 
sition, removing correlations on scales larger 
than half the box size. In addition, each box 
is assigned a bulk velocity along the line-of- 
sight according to the linear theory. 



5. We repeat steps 2-4 for all pairs of outputs 
in the descending order from the last output 
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at z — 4: {a 
simulation at 



0.2) until the start of the 
100. 



6. After arranging all the boxes, we calculate 
the CMB map on a square patch on the 
sky (corresponding to the image of our cu- 
bic computational box at 2: = 7.7, the epoch 
when about half the box is ionized) with a 
given resolution (we achieve the highest res- 
olution of 512^ pixels). 

This procedure allows us to construct an im- 
age on the sky. However, it is customary to 
present the CMB anisotropics on the sky in terms 
of the power spectrum Cg, which is the spherical- 
harmonic transform of the correlation function 
C{6). In the small-angle limit, the spherical- 
harmonic transform reduces to the Hankel trans- 
form of zeroth order, 



■I 



Ci = 2n / deejo{et)C{e), 



(12) 



where Jq is the zero order Bessel function. How- 
ever, because of the numerical errors and artificial 
ringing due to a finite box size, the power spec- 
trum calculated from equation (12) is noisy. An 
alternative method of calculating the power spec- 
trum is a direct Fourier transform of the CMB 
anisotropics on the sky, 



Q = Cax(|Arr), 



(13) 



where ^max is the angular size of a square patch 
on the sky (obviously, ^max l)j and 



is the Fourier transform of the CMB anisotropics. 

Figure 2 shows a comparison between the two 
methods for computing the power spectrum C^. 
One can see that the correlation function method 
produces an extremely noisy estimate of the power 
spectrum, whereas the one computed directly is 
well defined up to the maximal possible multipole 

Cax = \/2V = \/2Af-^, 

t'max 

where ^Ny is the Nyquist frequency of a square 
image on the sky and N is the number of pixels 
along one dimension in the image. 
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Fig. 2. — Comparison of the power spectrum of 
the CMB anisotropics computed by equation (12) 
{dotted line), and by equation (13) {solid line). 
In the latter case also shown are the statistical 
error-bars for the computed power spectrum. Two 
arrows show the angular size of the computational 
box and the angular resolution of our calculation 
respectively (for the 256^ pixelization). 
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Fig. 3. — Corrected {bold lines) and uncorrected 
{thin lines) power spectra at four different angular 
resolutions: 64^ pixelization {dotted lines), 128^ 
pixelization {short-dashed lines), 256^ pixelization 
{long-dashed lines), and 512^ pixelization {solid 
lines). Notice that our correction fully removes 
the effect of pixelization on small scales. 
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The process of pixelization - binning of the com- 
puted image on the square patch on the sky of a 
given size - produces small-scale smearing similar 
to the one produced by assigning density on a reg- 
ular mesh from a particle distribution in cosmo- 
logical A''-body simulations. However, in our case 
the precise form of this smoothing cannot be com- 
puted from the first principles, because the compu- 
tational box subtends somewhat different angular 
size at different redshifts. Instead, wc correct for 
the small-scale smearing empirically. Let Q be 
the true power spectrum. Since the smearing is 
caused by pixelization, the computed power spec- 
trum Ci will be a product of the true power spec- 
trum and the window function W which depends 
only on the ratio of the multipole £ and the maxi- 
mal possible multipole imax for a given pixel size, 
or, in other words, on the product of £ and the 
number of pixels along one dimension N. Thus, 

Ce = W{£/£^^)Ci, 

and if we take the ratio of the power spectra from 
two different pixclizations, this ratio will depend 
only on the window function, 



C2.e 



(14) 



Using equation (14), we find the following approx- 
imate form for the window function 



W{£/£„ 



10/3(^/Cax)']' 



1 + (^/Cax)^ 



(15) 



Figure 3 shows the uncorrected and corrected 
power spectra for four resolutions: N = 64, 128, 
256, and 512. One can see that the corrected 
power spectra agree with each other up to their 
respective ^max- We also notice that the power 
spectrum seems to continue as a power law with in- 
dex -3 to the smallest scales. Thus, we are unable 
to achieve convergence on small angular sc;alcs, 
which implies that even with resolution of 512^ 
pixels per 2.2 arcmin patch on the sky, there ex- 
ists structure on the unresolved scales. As we men- 
tioned above, this structure is due to small-scale 
density inhomogeneity in the ionized high density 
regions. Our simulation has a dynamical range 
of 4000, and thus in principle we can measure the 
power spectrum of the secondary anisotropies from 
our simulation up to Z 3 x lO'' (compared to 



I — 3.5 X 10^ for our highest 512^ resolution). How- 
ever, making CMB maps with resolution higher 
than 512^ is beyond the limit of computational re- 
sources available to us (doing a 4000^ image would 
take more computer time than the simulation it- 
self consumed). We thus conclude that the cut-off 
in the power spectrum at Z ^ 3 x 10"'' reported in 
Bruscoli et al. (2000) is not real and due to the 
finite angular resolution (256^) of their final CMB 
map. 

The outputs from our simulation are only 
stored at discreet moments in time (which we 
choose to parametrize with the cosmological scale 
factor a = 1/[1 -h z]). Integration over the path 
of a photon is done numerically using the out- 
puts from the simulation, and thus the numerical 
value of the integral depends on the interval Aa 
between the consequent outputs, converging to 
the exact value in the limit Aa 0. Of course, 
we cannot reach this limit as it would require an 
infinite number of outputs, and thus we must en- 
sure that our choice of Aa is sufficiently small to 
give an accurate value for the integral over the 
photon path. Figure 4 shows the power spectra 
of the CMB anisotropies for two choices of Aa. 
Both values give convergent results for the CMB 
anisotropies, and so we adopt Ao = 0.005 as our 
choice for this parameter (which is equivalent to 
having about 40 outputs until z = 4). It also 
becomes prohibitively expensive to calculate the 
integral with a smaller Aa. We also show in Fig. 
4 with the long-dashed line the CMB anisotropies 
calculated without randomization (i.e., in an un- 
realistic periodic universe), which are off by up to 
a factor of 10. 

We would like to note here that a recent similar 
analysis by Bruscoli et al. (2000) did not use ran- 
domization in calculating the CMB anisotropies, 
and by adopting a periodic universe, overesti- 
mated the correct result by a factor of 3-10. 

Finally, since our simulation is stopped at z = 
4, we need to estimate the contribution from the 
lower redshifts. Since the universe is fully ion- 
ized at z = 4, this contribution is entirely due to 
a homogeneous ionization fraction, Xg = 1 (the 
NL Ostriker-Vishniac effect). We can calculate it 
by extrapolating from higher redshifts. Figure 5a 
shows the rms CMB temperature anisotropy as a 
function of the final redshift of calculation z from 
our simulation (up to z = 4) with the solid line. 
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Fig. 4. — Power spectra of the CMB anisotropies 
pixelizcci on the 256^ image computed with two 
different values of the sampling frequency of nu- 
merical integration: Aa = 0.01 {dotted line), 
Aa = 0.005 {solid line). The dashed line shows 
a calculation with no randomization, i.e., with the 
periodic universe, which gives about a factor of 3- 
10 larger signal. Two arrows give the angular size 
of the computational box and the angular resolu- 
tion of our calculation respectively (for the 256^ 
pixelization). 
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Fig. 5. (a). The rms CMB temperature 
anisotropy as a function of the final redshift of 
integration. The solid line shows our simulation, 
which is stopped at z = 4, and the dotted line 
shows an extrapolation to z = assuming that 
perturbations do not change in comoving coor- 
dinates. (&). A comparison of the perturbation 
power spectrum as extrapolated to z = {solid 
line) and at 2: = 4 {dashed line). The difference 
between the two curves is representative of the un- 
certainty of our calculation. 
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One can sec that even aX z ~ A the rms tempera- 
ture anisotropy continues to rise (despite the fact 
that the universe is fully ionized by then), and 
thus we miss a portion of the total signal. In or- 
der to estimate the missing signal, we extrapolate 
the signal to ^; < 4 assuming that the density and 
velocity structure remains fixed in the comoving 
coordinates (in other words, taking the output of 
the simulation at 2; = 4 and adopting it for all 
lower redshifts). Then we find that 



AT 



(2 = 0) = 1.25 



AT 



RMS 



(2 = 4), 

RMS 



i.e., we miss about 20% of the signal. This number 
is likely to be an overestimate, as the gas pressure 
effects will erase the structure on spatial scales up 
to 500 kpc by 2 = 0, which corresponds to an- 
gular scales which are much larger than the ones 
dominating the signal. Nevertheless, since we can- 
not continue the simulation beyond 2: = 4, we use 
the extrapolated power spectrum as our final re- 
sult. The two spectra, the extrapolated one and 
the one directly computed from the simulation, are 
shown in Fig. 5b. The difference between the two 
is larger than, albeit comparable to, the statistical 
error due to the finite box size of the simulation, 
and is the largest uncertainty of our calculation 
on small angular scales (must less than the size of 
the box). On the angular scales comparable to the 
box size, an additional uncertainty is introduced 
by the missing large-scale power. 

3.3. Perturbations Outside of the Simula- 
tion Box 

Since the size of the simulation box is fixed to 
Ah~^ Mpc in comoving units, we are unable to cal- 
culate directly anisotropies due to perturbations 
on larger scales. However, we can include them 
analytically using the linear theory. Since we are 
not able to follow the ionization state of the gas on 
these scales, we assume that in the linear regime 
the ionization fraction is spatially homogeneous. 
This is not a bad approximation after all, since 
the size of a typical Hii region is smaller than the 
simulation box size, and thus on larger scales we 
average over a sufRciently large volume of the uni- 
verse. The latter implies that in the linear regime 
we can consider only the Ostriker-Vishniac effect 
- the NLOV effect is the same as the LOV effect. 



As we argue below, it is the dominant contribution 
to the total anisotropies on all angular scales. 

To calculate the pure LOV effect (linear fluc- 
tuations and homogeneous reionization), we use 
the approach described in Jaffe & Kamionkowski 
(1998). We use a small-angle Fourier-space ver- 
sion of Limber's equation which decouples the 
line-of-sight projection from the three-dimensional 
power spectra. As discussed above, we also assume 
the linear relationship between density and veloc- 
ity, uniform reionization, and Gaussian statistics. 
This reduces the calculation to a multiple integral 
over a product of the linear power spectrum of 
density perturbations with itself. 

To check our methods against this analytic cal- 
culation, we perform an additional "simulation" of 
the LOV effect in which we assume that the den- 
sity and velocity fields evolve only according to lin- 
ear theory, and we assign the ionization fraction at 
a given redshift uniformly over the computational 
box with the value given by the volume average 
ionization fraction in the full numerical simula- 
tion. We then make an image and calculate the 
CMB power spectrum in precisely the same fash- 
ion as in the full simulation. Thus, the spectrum 
of the anisotropies from a such "linear" simulation 
includes all the finite-box effects present in the full 
numerical simulation. 

The full nonlinear power spectrum is shown 
with the bold solid line in Fig. 6. In addition, 
we show the output from our linear "simulation" 
along with the analytic OV calculation, which 
match over all scales computed. We emphasize 
that this agreement requires the imposition of ap- 
propriate velocities on scales larger than the sim- 
ulation box, as discussed above. We also note 
that our full nonlinear prediction (bold solid line) 
matches the linear theory on the angular scale 
corresponding to the box size. This is of course 
is an artifact of our procedure (and the limited 
size of our simulation box), which does not in- 
clude nonlinear perturbations outside the simula- 
tion box. In reality, there exists additional power 
on scales which wc can not resolve, simply because 
the nonlinear scale at z = is about twice larger 
than our whole box. To account for this power, 
wc show the nonlinear power spectrum computed 
by Springel et al. (2000). They assumed precisely 
the same cosmological model and the redshift of 
reionization. In addition, their mass resolution al- 
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Fig. 6. — The spectrum of the CMB temperature 

anisotropy as calculated in linear theory [short- 
dashed line), from the full numerical simulation 
{hold solid line), and from the "linear" simula- 
tion {dotted line). Also shown with the long- 
dashed lines the nonlinear power spectra (for two 
box sizes) of kinetic SZ effect on large scales from 
Springel et al. (2000). 




Fig. 7. — Secondary temperature anisotropics in 
a 2.2' X 2.2' patch on the sky as generated dur- 
ing reionization {left panel) and with randomized 
phases (i.e., a realization of a Gaussian with the 
same power spectrum, right panel). The color 
bars show the color-scale correspondence, where 
the anisotropics are measured in jiK. 



most exactly corresponds to 1/8 of the total mass 
of our simulation box (corresponding to 1/2 the 
box size - the largest scale which can be resolved 
in our simulation). Thus, Springel et al. (2000) 
begin precisely at the scale where our calculation 
ends, and we will use their results to compliment 
ours below. (Conversely, this means that we un- 
fortunately cannot actually check the calculations 
against one another.) 

4. Results 

4.1. Maps and non-Gaussianity 

It is customary to characterize the CMB tem- 
perature anisotropics on the sky by the power 
spectrum Ci. The power spectrum is sufficient 
to fully describe properties of a physical quantity 
only if this quantity is Gaussian distributed. For 
the primary temperature anisotropics this is in- 
deed the case, since they are generated in the lin- 
ear regime. Reionization however is a non-linear 
process, and secondary anisotropics we consider in 
this paper are not necessarily Gaussian. 

This is illustrated in Figure 7, which shows 
the actual maps of the CMB anisotropics on the 
square patch on the sky in the simulation. The left 
panel shows the actual anisotropics on 2.2 arcmin 
patch on the sky (the angular size of our compu- 
tational box at z = 8, just before the full overlap 
of Hii regions). The right panel, however, shows 
anisotropics with the same power spectrum but 
distributed as a Gaussian. We notice that the ac- 
tual map is highly non-Gaussian; it is much more 
structured than the Gaussian version and contains 
higher fluctuations than would be expected in a 
Gaussian field. 

To illustrate this further, we show in Figure 
8 the one-point distribution of the actual tem- 
perature fluctuations, and a Gaussian distribution 
for comparison. We notice that the temperature 
anisotropics are non-Gaussian, and thus more in- 
formation than simply can be extracted from 
the actual distribution on the sky. In addition 
to temperature anisotropics on the sky, we also 
show in Fig. 8 the one-point distribution func- 
tion of spherical multipole amplitudes a^m, which 
are defined in the usual fashion by expansion over 
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spherical harmonics, 




Fig. 8. — One-point probability distribution of 
the temperature anisotropies {solid line), the one- 
point probability distribution of spherical multi- 
pole amplitudes (dashed line), and a Gaus- 
sian distribution {dotted line), all normalized to 
unit dispersion. 
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Fig. 9. — The electron density correlation func- 
tion Cee(?7, ?? + A?;, 9) as a function of angular size 
9 (measured in units of comoving distance Rg at 
this redshift) and A77 (also measured in units of 
comoving distance i?^) at three redshifts (values 
of 77): (a) z ~ 10, (6) z ^ 7, and (c) z = 5. Sub- 
sequently heavier shaded contours show increase 
in the correlation function by a factor of 10. The 
outmost contour corresponds to Cge = 1- 



— {9,<P)=J2<'^rnYr{9,cl,). 
i,m 

The one-point distribution of the amplitudes is 

actually quite similar to a Gaussian, which indi- 
cates that it is indeed the the phase correlations 
between different multipoles what drive the non- 
Gaussian patterns. In order to characterize this 
non-Gaussianity, one would have to measure the 
higher-order moments (bispectrum, trispectrum, 
etc.) at non-zero lags, or other quantities such as 
the Minkowski functionals, but that is outside the 
scope of this paper. 

4.2. Correlation Functions 

We now investigate the behavior of the correla- 
tion functions. In particular, the electron density 

correlation function Cee {Vi tV2,9) is of highest in- 
terest. We show in Figure 9 the electron density 
correlation function at three redshifts: before the 
overlap of Hii regions at z = 10, during the overlap 
at z = 7, and after the overlap at z — b. One can 
notice that there is no drastic change in the shape 
of the correlation function despite the qualitative 
change in the structure of the ionized regions; only 
the correlation length increases with time. This 
indicates that the signal is dominated by the ion- 
ized high density regions at all times, and is not 
sensitive to the ionization state of the low density 
intergalactic medium. We will elaborate on this 
below. 

Figure 10 shows both the electron density cor- 
relation function Cee and the velocity correlation 
function as a function of angle (measured 
as comoving distance at this redshift) at ^ = 7. 
We also show with the dotted line the quantity 
C^y/Cyy, which is much smaller than Cee if the 
cross-correlation between the electron density and 
velocity is not significant, and is equal to Cee is 
they are completely correlated. One can see that 
the cross-correlation only becomes significant on 
very small scales, about 0.1 arcscc {£ ^ 6 x 10®), 
and thus the usual assumption that velocity and 
electron density are uncorrelated is highly accu- 
rate on all scales of interest. 
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4.3. Kinetic vs Thermal SZ Effect 




Fig. 10. — The electron density correlation func- 
tion Ceei^, ^) {solid line) and the velocity corre- 
lation function Cyy{r],r],0) {dashed line) at z = 7. 
The dotted line shows the quantity C^^/Cyy. 
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Fig. 11. — Comparison of the power spectrum of 
the secondary CMB anisotropics in the computa- 
tional box generated by the kinetic SZ effect (eq. 
[1], solid line) and the thermal SZ effect (eq. [3], 
dashed line). 



It is interesting to compare the kinetic SZ effect 
given by equation (1) with the thermal SZ effect, 
given by equation (3). One can expect a priori 
that the kinetic SZ effect, generated in the gas 
moving with some 300 km/ s is about 1000 times 
more important (in terms of the power spectrum) 
than the thermal SZ effect, generated due to ther- 
mal motions in the lO'^ K hot gas. Indeed, this is 
demonstrated in Figure 11. Thus, for all practical 
purposes, the kinetic effect is the only one that 
needs to be taken into account on these scales. 
Of course, for clusters of galaxies, the situation is 
reversed and the temperature fluctuation induced 
by even the several- hundred km/s bulk motion of 
the cluster is dwarfed by the thermal effect, simply 
because the cluster gas temperature is some four 
orders of magnitude higher than the IGM temper- 
ature, whereas the bulk gas velocity in a cluster is 
only about a few times that of the IGM. We return 
to the comparison between the two effects below. 

4.4. The Ostriker-Vishniac Effect versus 
Patchy Reionization 

We now focus on a comparison between the NL 
Ostriker-Vishniac effect and patchy reionization. 
Figure 12 shows the comparison between the two 
correlation functions and , which are de- 
fined by equations (10). The dotted line, mark- 
ing the cross-correlation term in comparison with 
C^J^, only falls significantly below C^J^ at large 
radii {R > 0.1h~^ Mpc, where the inhomogeneous 
distribution of the ionized fraction due to expand- 
ing Hii regions make the NL Ostriker-Vishniac ef- 
fect and patchy reionization effect uncorrelated. 
However, as the thin dashed line shows, in the 
high density regions (i.e., on small scales) the ef- 
fect of patchy reionization is almost equivalent to 
having Xg = 1, i.e., to having almost all high den- 
sity regions ionized. This implies that the bulk 
of the signal comes from high density regions and 
correlations between them. 

In order to investigate the relationship between 
the two effects further, we have constructed a 
simulation which contains only the "Ostriker- 
Vishniac" effect by assigning a uniform ionization 
fraction (equal to the volume average at a given 
rcdshift) to the outputs of our full numerical sim- 
ulation, and performing line-of-sight integration 
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Fig. 12. — The electron density correlation func- 
tion CeeiViVi^^) for the patchy rcionization effect 
{C^J^, solid line) and for the NL Ostriker-Vishniac 
effect (C£^, bold dashed line) at ^ = 9 (at other 
redshifts the two correlation functions look quite 
similar except for the change of the characteris- 
tic scale). We also show with the dotted line the 

quantity {Cf^'^^f /C£^, which coincides with 
if cross-correlations are important. The thin 
dotted line shows the NL Ostriker-Vishniac cor- 
relation function rescaled by the quantity 

[(1 - Xe)/Xef. 



as described above. Thus, this "NLOV simula- 
tion" has the density and velocity structure of 
our full numerical simulation, and has the same 
evolution of the volume averaged ionization frac- 
tion, but has uniform ionization at all densities 
and spatial locations. In particular, the Patchy 
Rcionization correlation function is equal 

to zero in that simulation, and its NL Ostriker- 
Vishniac correlation function (C£^) is the same 
as in the full simulation. Thus, the NLOV simula- 
tion allows us to construct an image of the CMB 
anisotropics which correspond to just the 
correlation function, whereas the full simulations 
gives us anisotropics corresponding to the sum of 
the NL Ostriker-Vishniac and Patchy Rcionization 
effects. 

Figure 13 now shows the rms temperature 
anisotropy on the sky (for the 256^ pixelization) 
for the total effect, and separately for the NL 
Ostriker-Vishniac effect (as computed from the 
"NLOV simulation" ) and for the patchy rcioniza- 
tion effect (computed as the difference between 
the two). Note, that hy z = A almost all the sig- 
nal comes from the NL Ostriker-Vishniac effect, 
and the contribution of the patchy rcionization 
actually slightly declines at < 7, after the over- 
lap of Hii regions. This is due to the fact that 
after the overlap the topology of the ionized phase 
flips over: rcionization starts with isolated ion- 
ized regions expanding into the neutral medium, 
and ends with the isolated neutral regions being 
ionized from outside inward, which produces an 
anti-correlation between the ionized and neutral 
phases and leads to suppression of anisotropics 
due to patchy rcionization. 

Figure 14 presents now the main result of this 
paper: the power spectrum of the secondary CMB 
anisotropics from the kinetic SZ effect, from the 
thermal SZ effect, and from the patchy rcioniza- 
tion effect. For the former two, we combined our 
results with the calculations by Springel et al. 
(2000). In the range 3 x lO'' < ^ < 2 x lO'"^ our 
calculation for the kinetic SZ effect disagrees with 
the extrapolation from the Springel at al. (2000), 
and for the thermal SZ effect the range of dis- 
agreement is even larger. Since our results can 
suffer from the missing large-scale power at low 
redshift, and small scale extrapolation of Springel 
et al. (2000) is likely to be severely affected by lim- 
ited numerical resolution, it is impossible to give 
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Fig. 13. — The rms CMB temperature anisotropy 

as a function of the final rcdshift of integration for 
the NL Ostriker-Vishniac effect {dashed line), the 
patchy reionization effect {dotted line), and for the 
total effect (the sum of the two, solid line). 




Fig. 14. — The solid line shows the total ki- 
netic SZ power spectrum of the secondary CMB 
anisotropies, computed by combining our full non- 
linear results (the solid line from Fig. 6) for £ > 
10^ with the results of Springcl ct al. (2000) for 
10^ < £ < 3 X 10" and linear theory for £ < lO^), 
whereas the long-dashed line shows the same cal- 
culation for the thermal SZ effect. For both the ki- 
netic and thermal effects, the bold portions of the 
lines show the regions of the spectrum reliably cal- 
culated in this work (large £ tail) , in Springcl ct al. 
(2000), or in linear theory, and the thin line shows 
the regions of £ where existing predictions are un- 
reliable (and thus the curves we show are some- 
what arbitrary). Two dotted lines show the result 
for the patchy reionization, calculated from our 
simulation as the difference between the full sim- 
ulation and the "NLOV simulation" {bold dotted 
line), and an approximate prediction from Gruzi- 
nov & Hu (1998) {thin dotted line). Also shown 
is the spectrum of primary anisotropies for the as- 
sumed cosmological model {thin solid line). 
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a reliable prediction for the CMB power spectra 
in this range of angular scales, and we emphasize 
this fact by using the bold lines to show the calcu- 
lations where they are reliable, and by using the 
thin lines in the range of scales where consider cal- 
culations to be unreliable. 

We also show the prediction from an analyt- 
ical calculation of the patchy reionization effect 
by Gruzinov & Hu (1998) for our cosmological 
model. Note that this prediction is a factor of 
10^ lower than the one given in Gruzinov & Hu 
(1998) for the standard CDM model, because they 
adopted the value for the rms velocity at z = of 
1200 km/ s, whereas in our model this number is 
about 3 times smaller, in agreement with obser- 
vations (Baker, Davis, & Lin 2000). There is an 
additional reduction in the amplitude of the effect 
by a factor of 10 due to the redshift of reionization 
being 7 in our simulation compared to an adopted 
value of 30 in Gruzinov & Hu (1998). Including 
all the relevant factors, we notice that Gruzinov 
& Hu (1998) estimate gives roughly the right pre- 
diction for the patchy reionization effect on large 
scales. Of course, since the analytical estimate of 
Gruzinov & Hu (1998) does not include the den- 
sity structure on small scales, they are unable to 
reproduce the anisotropies on small angular scales. 

4.5. Dependence on the Redshift of Reion- 
ization 

Up to know we assumed a fixed value for the 
redshift of reionization, since we only used one 
simulation. In order to investigate the dependence 
of the power spectrum on the redshift of reioniza- 
tion properly, we would need to run several sim- 
ulation with different reionization histories, which 
is not realistic at the current moment due to re- 
quired computational resources. However, in or- 
der to illustrate the dependence of anisotropies on 
the redshift of reionization, and using the fact that 
at zrei = 7 the NL Ostriker-Vishniac dominates 
the contribution to the anisotropies due to Patchy 
Reionization, we have computed three additional 
"NLOV" simulations where we changed the red- 
shift of reionization by simply rescaling the vol- 
ume averaged ionization fraction as a function of 
redshift. In other words, if ne.o(a) is the electron 
density and Xe{a) is the volume averaged ioniza- 
tion fraction as a function of the scale factor for the 
original "NLOV" simulation, for rescaled "NLOV" 



simulations we assumed 

Xe(8aREia) , . 
ne[a) = — — - — ne,o(a), 
Xe{a) 

where the factor 8 is simply 1/0.125 = 1/aREi = 
1 + 2r,ei for our simulation. 

Figure 15 shows the CMB power spectra for 
four different values for the redshift of reioniza- 
tion in our rescaled "NLOV" simulations. We note 
that in linear theory, 5 a and v ^ a^/^ at high 
redshift, and therefore the total anisotropy 
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ribdt ~ log(t) 



goes only logarithmically as a function of the mo- 
ment of reionization (which is the upper limit in 
the integral [1] for a model with sudden reioniza- 
tion). The dependence that we find in the non- 
linear calculation is somewhat steeper than that 
from the linear theory, but not by much. 

We also note here that the analytical model of 
patchy reionization by Gruzinov & Hu (1998) pre- 
dicts that the spectrum of the CMB anisotropies 
increases as (1 -|- ^;rei)^/^-RjjiI' where -Rjjii 

characteristic comoving size of Hii regions, which 
is likely to go down as the redshift of reioniza- 
tion increases. Thus, the total dependence on 
the redshift of reionization of the patchy reioniza- 
tion effect is somewhat slower than (1 + zrei)^^^ 
but perhaps not as slow as the nearly logarith- 
mic dependence of the NL Ostriker-Vishniac ef- 
fect. Thus, if reionization occurred early, the role 
of patchy reionization (as compared to the NL 
Ostriker-Vishniac effect) would be greater, albeit 
the current observational limits of zrei ^ 30 will 
not allow a regime when patchy reionization be- 
comes dominant. 

5. Conclusions 

We use a cosmological numerical simulation 

that models time-dependent and spatially-inhomo- 
geneous cosmological reionization to compute sec- 
ondary CMB temperature anisotropies in a rep- 
resentative CDM+A cosmology. We compare our 
numerical results with analytical calculations and 
compliment them with the numerical simulations 
of Springel et al. (2000) (whose mass resolution al- 
most exactly compliments ours, and accidentally 
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Fig. 15. — The power spectra of the secondary 
CMB anisotropics computed in the "NLOV" sim- 
ulation with four different values for the redshift 
of reionization: zrei = 5 {dotted line), zrei = 7 
{short-dashed line), z^^i = H {long-dashed line), 
and z^Ei = 19 {solid line). 



both simulations have the same values of cosmo- 
logical parameters). We thus are able to compute 
the power spectrum of secondary anisotropies for 
the kinetic SZ effect and for the patchy reioniza- 
tion effect over a wide range of angular scales. 
The thermal SZ effect, dominated at medium £ 
by massive clusters of galaxies, is more difficult 
to compute on all scales due to the extreme non- 
linearity of the problem. 

We find that the role of patchy reionization 
per se (i.e., evolution of Hii regions around the 
sources of ionization) in generating the secondary 
anisotropies is subdominant. This conclusion is 
however only valid when the sources of reioniza- 
tion are clustered on scales comparable to the mat- 
ter correlation length, in which case one can ex- 
pect to find a source (or multiple sources) in al- 
most every dark matter halo. If, in the oppo- 
site case, the sources of reionization arc rare (like 
bright quasars), the relative role of patchy reion- 
ization will be greater, as many high density re- 
gions will remain neutral for a long time, and may 
not be negligible. The latter scenario is, however, 
less likely given the current observational limits 
on the abundances of bright quasars (Fan et al. 
2000b), and thus CMB observations on small an- 
gular scales may not be as useful as hoped to study 
the details of reionization. 

A real window onto reionization is provided by 
the detailed morphology of the fluctuations, that 
is, the non-Gaussianity discussed above. The an- 
gular scale of the fluctuations traces the three- 
dimensional separation scale of the ionized regions. 
Of course, a full understanding of this will re- 
quire not just a detection of broad-band fluctu- 
ation power at very small angular scales, but a 
detailed, high signal-to-noisc map of these very 
small-scale fluctuations with an angular resolution 
of arc seconds, which is definitely some years away. 

While the details of reionization may not be 

easily discernible from the broad-band power spec- 
trum, the amplitude of the fluctuations, which 
is directly related to the redshift of reioniza- 
tion, surely is. However, as our comparison with 
Springel et al. (2000) shows, on medium angu- 
lar scales (/ < 10^) the signal is dominated by 
lower redshifts (z ^ 2 — 4). and only on small 
angular scales (several arc seconds) the amplitude 
of the anisotropies is related to the redshift of 
reionization. In addition, this dependence is quite 
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weak, unless rcionization occurred close to the 
currently allowed observational limit of zrei ^ 30. 
Nevertheless, observations on those angular scales 
can put significant constraints on the rcdshift of 
reionization. Such arcsecond-scale CMB observa- 
tions are best accomplished with interferometers, 
and indeed the recent development of compact 
submillimeter-wave arrays has paved the way for 
such observations. For example, the ATCA tele- 
scope finds a limit of {e{e+l)Ci/{2TT)) < 2 x IQ-i" 
at f ~ 4000 (Subrahmanyan et al. 2000), well 
above the values presented here. In addition, 
these observations are inevitably plagued by point 
source confusion, especially dangerous when the 
signal also has a compact galaxy-like component. 
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